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1.  INTRODUCTION 


Since  the  creation  of  a  light  weight  composite  for  the  space  industry,  com¬ 
posite  materials  have  now  been  widely  used  in  many  areas  of  the  industry.  Due 
to  the  geometrical  and  material  discontinuity  in  the  composite,  stress  singular¬ 
ities  often  occur  resulting  in  the  failure  of  the  composite.  In  order  to  under¬ 
stand  the  failure  mechanism  and  to  help  designing  a  better  composite,  it  is  es¬ 
sential  to  analyse  the  nature  of  stress  distribution  near  a  singular  point  in 
the  composite.  One  of  the  reasons  for  the  slow  progress  in  designing  a  more  du¬ 
rable  composite  is  the  lack  of  rigorous  stress  analyses  near  a  stress  singulari¬ 
ty,  especially  when  the  singularity  is  a  three-dimensional  one.  Even  though 
many  investigations  have  been  carried  out  on  the  nature  of  stress 
singular i ties  in  2-dimensional  elasticity  problems,  only  a  relatively  little 
work  has  been  done  in  the  area  of  stress  singularities  in  general  3~D  problems. 
This  may  be  partly  due  to  the  complexities  of  the  problem  itself.  Nevertheless, 
in  practice,  both  in  the  areas  of  fracture  mechanics  and  general  stress  analysis 
a  knowledge  of  the  stress  singularities  would  be  quite  useful.  In  an  excel¬ 
lent  review  article  on  the  three  dimensional  stress  problem  for  a  cracked  plate 
Sih  [1]  discusses  at  length  the  difficulties  associated  with  the  task  of  obtain¬ 
ing  an  analytical  solution.  Most  of  the  time  in  analytical  studies  special  ma¬ 
terial  properties  and  geometries  that  make  it  possible  to  reduce  the  problem  to 
two  dimensions  in  mathematical  terms  have  been  employed.  An  extensive  collec¬ 
tion  of  such  analysis  is  presented  in  [2]  and  further  discussions  on  3~D  analy¬ 
ses  can  be  found  in  [3].  More  recently  Parihar  and  Keer  [I*, 5]  have  presented 
analytical  solutions  for  wedge-shaped  cracks,  inclusions  and  stamps  in  isotropic 
materials.  Such  special  problems  serve  as  extremely  useful  'bench  mark'  tests 


1 


for  other  more  genera)  methods  such  as  the  one  that  would  be  described  herein. 
However,  in  the  more  general  case  of  a  three  dimensional  singularity  in  a  com¬ 
pletely  anisotropic  material,  obtaining  an  analytical  solution  seems  quite  dif¬ 
ficult  except  for  special  cases  [6]  and  numerical  methods  have  to  be  employed. 
The  present  study  is  aimed  at  proposing  one  such  numerical  method,  viz.  a  finite 
element  scheme. 

A  general  numerical  procedure  for  determining  3“D  singularities  was  first  de¬ 
veloped  by  Bazant  [7]  for  harmonic  functions  and  by  Bazant  and  Estenssoro  [8,9] 
for  isotropic  elastic  materials  who  used  it  to  study  the  singularities  at  the 
intersection  of  free  surface  and  cracks.  Their  results  for  a  crack  orthogonal 
to  a  free  surface  were  in  good  agreement  with  the  results  obtained  via  a  semi- 

analytical  method  by  Benthem  [10,11],  However  no  solutions  are  available  for 

« 

anisotropic  materials. 

The  present  finite  element  scheme  is  an  extension  of  that  developed  by  Bazant 
and  Estenssoro  [9]  to  incorporate  general  anisotropic  elastic  materials.  In 
Section  2  we  present  the  general  formulation  of  the  finite  element  scheme. 
Testing  of  the  method  against  some  known  results  will  be  described  in  Section  3> 
In  Section  4  it  is  applied  to  study  a  problem  in  laminated  composite  materials 
and  finally  in  Section  5  we  use  it  to  study  the  possible  occurrence  of  stress 
singularities  in  other  3"0  corners. 

Let  (r,6,<t>)  refer  to  spherical  co-ordinates  and  consider  the  region  near  the 
apex  of  a  conical  wedge  (a  notch  or  a  rigid  inclusion)  in  which  the  apex  is  the 
origin  of  the  co-ordinates.  Fig.  1.  The  lateral  surface  S2  of  the  conical 
wedge  is  assumed  to  be  generated  by  radial  lines.  The  cross  section  of  the  con¬ 
ical  wedge  at  any  constant  radius  is  denoted  by  Sf.  The  boundary  of  $  , 


which  specifics  the  shape  of  Sjf  is  given  by  the  equation  \ ■  0.  Depend¬ 
ing  on  whether  a  notch  or  a  rigid  inclusion  is  being  analysed  the  homogeneous 
boundary  conditions  applicable  on  S?  will  be  either  traction  free  condition  or 
rigidly  clamped  condition  respectively.  Since  we  are  interested  in  the  possible 
stress  singularities  near  the  apex,  the  boundary  conditions  on  S,  need  not  be 
specified. 

For  the  purpose  of  an  asymptotic  analysis  for  possible  stress  singularities 
at  the  apex  of  the  conical  region  we  seek  non-trivial  displacement  solutions  of 
the  form 

u  ■  r^  f  |  (1) 

v  ■  rX  g  / 

w  ■  rx  h  J 

where  (u,v,w)  are  the  components  of  displacement  in  ( r,6,<t> )  directions  respec¬ 
tively.  This  displacement  field  is  required  to  satisfy  equations  of  equilibrium 
within  the  conical  region  and  the  relevant  homogeneous  boundary  conditions  on 
the  lateral  surface  S2. 

A  value  of  X  that  would  satisfy  these  conditions  is  called  the  eigenvalue  and 
the  corresponding  functions  g(0,<p)  and  h(0,<£)  are  called  the  eigenfunc¬ 
tions.  There  are  infinitely  many  such  eigenstates  [12].  But  in  the  present 

case  our  attention  will  be  confined  to  those  eigenstates  that  would  lead  to 
stress  singularities  at  the  apex.  When  displacements  are  of  the  form  given  by 
equation  (1)  the  stresses  are  proportional  to  r^-1.  Hence  a  stress  singulari¬ 
ty  occures  at  the  origin  when  X  <  1.  On  the  other  hand,  for  the  strain  energy 
to  be  bounded  at  the  origin  we  require  X  >  -J5.  However,  X  <  0  implies  that  the 
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displacements  are  unbounded  at  r  ■  0  which  is  physically  unrealistic  unless  a 
concentrated  load  is  applied  at  the  wedge  apex.  Therefore  we  are  primarily  in¬ 
terested  in  eigenvalues  X  in  the  range  0  <  X  <  1.  The  most  important  would  be 
the  smallest  X  (i.e.,  closest  to  0)  which  would  give  rise  to  the  strongest  sin¬ 
gularity.  If  X  is  complex,  we  are  interested  in  the  case  where  0  <  Re  (X)  <  1. 
During  this  procedure  no  particular  boundary  conditions  are  prescribed  at  the 
surface  S,  of  the  cone.  The  tacit  assumption  is  that  by  superposing  a  suffi¬ 
cient  number  of  eigenstates  they  can  be  satisfied. 

In  order  to  evaluate  X  (  and  f  (6,<t>),  g  and  h(0,4>)  )  that  would  satisfy 

the  above  conditions,  Bazant  and  Estenssoro  [9]  developed  a  finite  element 
scheme  based  on  a  variational  principle  involving  those  quantities.  Even  though 
they  considered  only  isotropic  materials,  with  a  few  changes  in  the  details 
their  approach  can  be  used  for  general  anisotropic  materials.  In  the  next  sec¬ 
tion  we  shall  present  the  general  formulation.  In  [9]  it  has  been  shown  that 
the  variational  principle  does  not  yield  a  minimum  principle  and  that  the  re¬ 
sulting  system  of  equations  for  the  determination  of  X  is  non  symmetric.  Those 
proofs  are  valid  for  the  present  formulation  also. 
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2.  FORMULATION 


2.1  VARIATIONAL  PRINCIPLE 

Normally  problems  of  elasticity  can  be  solved  by  the  use  of  stationary  energy 
principles.  In  the  present  case  too  we  can  base  our  solution  procedure  on  them. 
However,  some  special  considerations  and  procedures  are  required  due  to  the  fact 
thit  the  solutions  we  seek  are  not  expected  to  satisfy  al 1  the  boundary  condi¬ 
tions  of  the  problem;  they  need  satisfy  only  the  near  field  boundary  conditions. 

For  the  conical  solid  shown  in  Fig.  1  the  two  boundaries  are  r  *  r0  and 
i P(6,<t>)  -  0  and  are  denoted  by  St  and  S2f  respectively.  Our  intention  is  to 

find  a  set  of  eigenfunctions  for  the  displacement  field  in  the  form  of  equation 
(1)  that  would  satisfy  the  necessary  boundary  conditions  near  the  origin  -  i .e. 
on  S2.  Nothing  is  said  about  the  boundary  conditions  at  the  far  end  -  i .e.  on 
$,.  The  hope  is  that  by  superposing  a  sufficient  number  of  eigenfunctions, 
boundary  conditions  on  S,  may  be  satisfied.  In  fact  the  cone  shown  in  Fig.  1 
should  be  considered  only  as  one  part  of  a  boundary  value  problem  and  matching 
the  solution  for  the  cone  with  that  for  the  adjoining  region  on  S,  is  not  of 
immediate  concern. 

The  eigenfunctions  we  seek  are  therefore  required  to  satisfy  the  governing 
differential  equations  and  only  some  of  the  boundary  conditions  of  the  complete 
elasticity  problem  for  the  cone.  It  is  due  to  this  reason  that  some  modifica¬ 
tions  to  the  stationary  energy  principles  are  needed  for  the  present  problem. 
These  modifications  were  first  proposed  by  Bazant  and  Estenssoro  [9].  We  too 
follow  their  procedure  with  a  slightly  different  presentation  of  the  rationale 
behind  the  modification. 

Once  again  consider  the  conical  region  shown  in  Fig.  1.  When  it  is  subjected 
to  some  displacement  field  a  certain  traction  will  result  on  the  boundary  S  . 
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Let  that  traction  be  denoted  by  p  a  (pr,pg,p^)  . 


If  one  were  to  solve  the 


I 


complete  problem  of  which  this  cone  is  only  a  part,  one  would  have  to  match 
these  tractions  and  displacements  on  S,  with  those  of  the  adjoining  region. 
But  for  the  purpose  of  our  local  analysis  consider  instead  the  cone  to  be  sub¬ 
jected  to  external  traction  p  on  S(.  Then  the  cone  should  be  in  equilibrium 
and  the  displacement  field  should  satisfy  the  equations  of  equilibrium  inside 
the  cone  and  the  relevant  boundary  conditions  on  S,  and  S2.  For  this  situ¬ 
ation  we  can  write  the  principle  of  minimum  potential  energy  as 


<$U  -  /  (pr<5u  +  P0<5v  +  p^iSw)  r*sin  8  d0dd>  *  0  (2) 

S1 

where 

U  m  f  <p  drd0d<£  (3) 

V 

is  the  total  strain  energy  in  the  volume  V  of  the  cone  and 


$  *  E  r’sin  fl  (1») 

where  E  is  the  strain  energy  density.  o  indicates  the  variations  due  to  kine¬ 
matically  admissible  variations  in  the  displacement  field  (u,v,w) .  E  is  a  func¬ 
tion  of  the  strains  «jj  and  they  in  turn  are  functions  of  the  displacements 
and  their  first  derivatives.  Therefore  we  get 


<5U  *  /  (  'fc.Su  +  5ur  +  . +  <t>w  5wA)  drd0d<i>  (5) 

V  r  <t>  w 

where  the  subscripts  r,9,4>  denote  partial  differentiation.  4>u  refers  to  the 

partial  derivative  of  $with  respect  to  u,  assuming  u,  ur,  uq, . w^  to 

be  independent  variables  and  so  on. 

Equation  (2)  implies  all  the  equations  associated  with  the  problem  of  elas¬ 
ticity.  But  the  eigenfunctions  need  not  satisfy  the  far  field  boundary  condi- 


tions  on  S i .  Therefore  we  will  modify  (2)  into  a  form  that  would  allow  those 
boundary  conditions  to  be  separated  out  from  the  others.  Substitute  (5)  into 
(2)  and  integrate  the  terms  containing  <5ur,  <5vr  and  <5wr  by  parts  in  the  r 
direction  to  yield: 


oU  +  f  [  $  5u  +  4>  5v  +  $  ow  1  d6d<t> 

S,  L  ur  vr  wr  J 

"  (  Pr$u  +  PgSv  +  p^Sw  )  rJsin  6  d6d<t>  - 


where. 


—  b 

5U  -  /  [  {*  -  — —  (<i>  )}  Su  +  <f>  <5ufl  +  *  Suz 

V  i  u  dr  ur  -  u<t>  * 

+  {%-  —  5v  +  *V05v0  +  S^Sv^ 

+  {V  i7(V}  5w  +  <Vws  +  \^<t>  ]  drdfld* 


To  evaluate  $Vr  and  4>W[_  first  evaluate  EU(_,  EVf. 

and  EWr<  Note  that  3E/d«jj  -  a jj  where  Ojj  is  the  stress.  Since 
«jj  is  related  to  u,  ur>...  through  the  strain  -  displacement  relations  in 
spherical  co-ordinates,  Ey^,...  can  be  easily  evaluated  using  the  chain 
rule.  We  then  obtain  » • • .  by  the  use  of  (4).  When  these  algebraic  ma¬ 
nipulations  are  carried  out  we  get 


*ur  -  r2  arr  sin  6 
%  m  r2  arQ  sin  6 
*wr  ■  r2  ffr<<>  *in  e 


Substitution  of  (8)  into  (6)  yields 
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6U  +  /  [  (<7rr-pr)  Su  +  (ffpg-pg)  Sv  +  (a^-p^)  6w  ]  r  sin  9  d66<t>  -  0 

S  i 


which  implies 


5U  -  0  , 


ffrr  "  Pr 
CTr 9  “  Pe 
ar<t>  *  Ptf 


Equation  (11)  is  in  fact  the  traction  boundary  condition  on  S,.  Since  we 
started  with  a  variational  statement  that  ensured  the  satisfaction  of  all  the 
boundary  conditions  (i.e.  on  both  $1  and  S2)  and  the  equations  of  equilibri¬ 
um,  we  conclude  that  (10)  is  a  variational  statement  that  would  ensure  the  sat¬ 
isfaction  of  the  equations  of  equilibrium  in  V  and  the  boundary  conditions  on 
S2  only.  Therefore  the  eigenfunctions  we  seek  should  satisfy  the  variational 
principle  (10).  It  should  also  be  noted  that  in  this  derivation  no  particular 
assumption  was  made  regarding  the  constitutive  relation  for  the  elastic  material 
except  that  it  possesed  a  strain  energy  density  function. 

Let  the  stresses  and  the  strains  be  represented  by  1-D  arrays  as  follows: 


ai  "  ffrr •  aj  *  aed>  “  a44> 
ff4  "  °6<t>>  *5  *  at  "  ar9 


•^  ‘  €rr • 


%  *  2  «5  "  2  ‘r*;  «,  ■  2  •re 

The  material  constitutive  relation  is 


“  Cij«j 


where  C ;  is  the  material  stiffness  which  satisfies 


In  equation  (14)  and  in  the  following  sections  repeated  indices  imply  summation 
unless  otherwise  stated. 

As  explained  earlier,  4>u,  4>u  ,  .  .  .  are  evaluated  by  using  equation  (4) 

and  treating  E  as  a  function  of  «jj  and  then  using  the  chain  rule.  The  re¬ 
sults  are: 


<J> 


$ 

*  i  i 


u0 


<j> 


J</> 


$ 


* 


* 


v0 


* 


4* 


<f> 


we 


4> 


W, 


<t> 


r  (o2  +  a3)  sin  0 
r2  o1  sin  8 
r  a#  sin  8 

r  a5 

r  (a3  cos  6  -  cr#  sin  Q) 
r2  at  sin  8 
r  «t2  sin  8 

r  ^4 

-r  (<74  cos  8  -  o5  s  i  n  9) 
r2  <75  sin  8 
r  a A  sin  8 


\ 


J 


(16) 


Now  we  can  use  the  expressions  in  equation  (1) 
displacement  relations  yield: 

«  j  *  r^- 1  t j  ,  i *1 , . . , 6 


for  u,  v  and  w.  The  strain 


(17) 


where 


*  a4  sin  9  I 

J 

Equation  (21),  when  substituted  into  (10), 
pie  for  the  present  problem  as 


(23] 
cont  ’d 


yields  the  final  variational  princi- 


5U  -  0  for  all  5f,  <5g,  Sh  (2it) 


Bazant  and  Estenssoro  [8,9]  who  investigated  the  3"D  stress  singularities  in 
isotropic  elastic  materials  using  this  approach  have  given  the  expressions  for 
F,  F17,...  for  isotropic  materials.  We  have  checked  and  verified  that  when 
Cjj  is  specialized  to  isotropic  materials,  the  expressions  given  in  (23)  do  in 
fact  reduce  to  those  given  in  [8,9]  (after  minor  printing  mistakes  are  corrected 
[9]). 

2.2  DESCRETIZEO  FORMULATION, 

To  derive  the  expressions  for  a  descretized  (finite  element)  formulation  let 
the  area  containing  material  bounded  by  the  curve  -  0  on  the  (6 ,4>)  plane 

be  denoted  by  A  and  let  it  be  sub-divided  into  an  n  number  of  finite  elements. 
The  values  of  f  (0,d>)  ,  g  and  h(0,<£)  are  taken  as  the  nodal  degrees  of  free- 

dom.  We  can  express  the  integral  5U  of  equation  (22)  as  the  sum  of  the  same  in¬ 
tegral  taken  over  each  element  (provided  that  the  continuity  conditions  dis¬ 
cussed  below  are  satisfied).  Hence 

5U  -  2  (25) 

where  <5U ^  refers  to  the  result  of  integrating  the  same  integrand  as  that 
given  for  <5U  in  equation  (22)  over  the  area  of  the  mth  element,  A  (m) .  2  in 

the  above  equation  as  well  as  in  the  following  steps  implies  summation  over  all 


1 1 


n  elements.  Examination  of  the  terms  in  the  integrand  of  equation  (22)  shows 
that  for  equation  (25)  to  hold  f,  g  and  h  should  be  C°  over  A  (i.e.  only  the 
functions  themselves  need  be  continuous,  not  the  derivatives).  This  justifies 
the  choice  of  the  values  of  the  functions  as  the  nodal  degrees  of  freedom.  The 
intra-element  interpolations  also  need  satisfy  only  function  continuity  along 
the  inter-element  boundaries. 

Let  £  be  the  vector  of  nodal  degrees  of  freedom  associated  with  a  typical  el¬ 
ement  and  define  interpolation  functions  L(0,</>),  M  (6,<t>)  and  N(0,</>)  by 


f  (8,*)  -  L'  (0.*)  X 
g(0,d>)  -  M1  (0,*)  X 
h  (0  ,</>)  -  )  X 


In  equation  (26)  as  well  as  in  the  following  steps  superscript  i  on  L,  M  or  N 
and  subscript  i  on  X  refer  to  the  i**1  entry  in  each  vector.  As  has  already 
been  stated  repeated  indices  implie  summation. 

Substitution  of  equation  (26)  into  (18)  gives 

«i  "  EjjXj  (27) 


with  the  2-0  array  E j j  defined  as 


E2j  “mJ0.+  L\ 

Ejj  ■  lJ  +  MJcot  0  +  /sin  0  l 

E4j-  -  -  NJ  cot  0  +  /sin  0  | 

E5j  -  (X— 1 )  NJ  +  /sin  0  I 

E6j  -  (X-l)  MJ  +  lJe  ' 

where  subscripts  0,  <t>  on  L-i ,  M-i  and  NJ  refer  to  differentiation  with  re¬ 
spect  to  that  variable.  When  equation  (27)  is  substituted  into  (20)  we  get 


where  the  2-D  array  Sjj  is  given  by 
Sij  "  cikEkj 

Now  by  substituting  equation  (29)  into  (23)  we  can  define 


wi  th 


pi 

?6\ 
P<i>  j 

Qi 

Q0i 

Q^i 

R; 

R0i 

R^i 


C  S2 i  +  S^j  -  (X+DSj,  3  sin  9 
S6i  sin  9 

S5i 

Sjj  cos  9  -  (X+2)S$j  sin  9 

S2 f  sin  9 

Hi 

-  S^j  cos  9  -  (X+2)Sj;  sin  6 

S4 j  sin  9 

s3i 


(30 


(31 


Using  equations  (26)  and  (31)  in  (22)  and  integrating  over  the  area  of  the  ele 
ment  A  we  finally  get 


X? 


where 

Kij)  "  ft(  )  {  Pjli  +  Pj$L^  +  ^  +  QjMi  +  +  ^  (34) 

+  RjN '  +  R$N£  +  RjN^  }  6d6<t> 

In  matrix  notation  (33)  is 

5U(m)  -  5XtK^X  (35) 

Let  X*  be  the  vector  of  global  degrees  of  freedom  and  let  X  be  related  to 

'V  rs* 

*  through  the  connectivity  matrix  B M  so  that 

X  -  B(m)X*  (36) 

When  this  is  substituted  into  equation  (35)  we  get 


5U  (m)  »  6X*TB  ("*)  \  (ffl)  B  ^  X* 

-V  A.  /V  <V 

(37) 

With  equations  (25)  and  (37).  the  governing 

expressed  as 

variational  principle  (24) 

can  be 

<5U  -  2  <5U  (m)  -  5X*TK  X*  -  0  for  all 

5X* 

(38) 

where 

*** 

K  -  2  B  (m)  TK  M  B  M 

^  ^  <v  *w 

(39) 

For  equation  (38)  to  hold  for  all  5X*  we  require 

K  X*  -  0  (1*0) 

/w  (V  <V 


Equation  (40)  identifies  K  as  being  similar  to  the  'global  stiffness  matrix'  in 
standard  finite  element  analysis  and  it  follows  that  K^m)  is  similar  to  the 
'element  stiffness  matrix'.  Hereafter  K  and  K  (m)  will  be  referred  to  by  those 


names . 


It  should  be  noted  that  equation  (39)  is  only  symbolic  and  in  practice 
the  formation  of  K  from  j<(m)  is  acompl i shed  through  the  usual  assembly  proce¬ 
dures. 

Elements  of  K  depend  on  X  and  to  derive  a  characteristic  equation  for  X  we 
note  that  for  equation  (40)  to  have  a  non-trivial  solution 

|  K  |  -  0  (41) 

This  provides  the  e i genequat ion  for  X. 

2.3  JHE  ELEMENT 

The  expressions  needed  for  evaluating  the  element  stiffness  matrix  have  been 
given  in  the  previous  section.  For  actual  numerical  implementation  of  the 
scheme  one  has  to  decide  on  a  particular  2-D  element.  Since  only  the  functions 
f,  g  and  h  are  required  to  be  continuous,  any  element  type  normally  used  in  2-D 
stress  analysis  can  be  picked  for  this  purpose. 

In  [8,9]  a  4-node  quadrilateral  based  on  the  same  variational  principle  has 
been  used  to  investigate  3“0  stress  singularities  in  isotropic  materials.  The 
results  showed  fairly  large  errors  and  meshes  with  a  substantial  number  of  ele¬ 
ments  were  required  to  achieve  a  reasonable  degree  of  accuracy.  However  it  was 
possible  to. exploit  the  convergence  pattern  of  the  results  to  extrapolate  them 
and  find  accurate  solutions.  When  a  symmetry  of  the  geometry  and  of  the  materi¬ 
al  property  exist,  one  can  make  use  of  the  symmetry  by  considering  only  one  half 
of  the  problem.  But  in  the  case  of  a  general  anisotropic  material  such  symmetry 
would  be  lacking  and  the  full  problem  will  have  to  be  analysed.  This  will  make 
it  difficult  to  use  very  fine  meshes  and  therefore  it  was  decided  to  use  a  high¬ 
er  order  element  for  the  present  study.  It  would  enable  reasonable  accuracies 
to  be  realized  without  having  to  use  a  large  number  of  elements. 


The  chosen  element  is  an  8-node  isoparametric  quadrilateral  with  nodes  at  the 
four  corners  and  the  mid-sides.  The  interpolations  are  quadratic  and  curved 
sides  are  admissible.  Following  standard  2-0  finite  element  procedures  all  in¬ 
tegrations  are  performed  numerically  on  a  parent  plane  onto  which  the  element  is 
mapped  as  a  unit  square.  This  element  is  now  standard  and  details  of  the  inter¬ 
polation  functions  and  numerical  integration  are  found  in  any  standard  text  book 
on  finite  elements,  (for  example,  see  chapters  7  and  8  of  [14]). 

Numerical  integration  is  performed  using  Gauss  quadrature  rule  on  a  grid  of 
ng  x  ng  integration  stations.  In  a  general  anisotropic  material  the  ele¬ 
ments  of  the  material  stiffness  matrix  Cjj  also  are  functions  of  9  and  <t>  and 
at  every  integration  station  they  too  have  to  be  evaluated.  The  transformation 
is  that  correspond i ng  to  a  4th  order  tensor  and  the  transformat i on  matrix  con¬ 
tains  terms  like  cos  6,  cos  d>,  sin  6  and  sin  d>.  Therefore  Cjj  will  display  a 
strong  dependence  on  these  terms  and  their  higher  powers.  The  terms  Ejj  also 
contain  these  trigonometric  functions.  Examination  of  the  expressions  involved 
show  that  terms  like  1/sin  6  also  are  present.  This  leads  to  some  difficulties 
numerically  in  evaluating  the  integrals  accurately  especially  in  the  region  near 
the  pole  of  the  coordinate  system  where  9  -  o.  Since  a  prior  estimate  of  the 
number  of  integration  stations  is  not  possible,  the  choice  of  integration  sta¬ 
tions  has  to  be  based  on  numerical  testing.  In  the  following  sections  where 
test  and  example  problems  are  presented,  we  will  illustrate  the  effect  of  vary¬ 
ing  the  number  of  Gauss  integration  stations  and  discuss  the  choice  of  optimal 
number  of  integration  stations. 


The  domain  occu- 


2.4  EVALUATION  OF  X 

The  details  of  computing  X  for  a  given  problem  is  simple, 
pied  by  the  material  in  the  {9,4>)  plane  is  subdivided  into  a  finite  element 
mesh.  For  a  particular  value  of  X,  element  stiffnesses  are  computed  and  are 
then  assembled  into  the  global  stiffness  matrix  K.  Boundary  conditions  are  im¬ 
posed  by  eliminating  rows  and  columns  corresponding  to  degrees  of  freedom  pre¬ 
scribed  to  be  zero  in  the  global  stiffness  matrix.  Finally  |  K  |  is  computed. 
The  eigenvalue  X  (for  which  )  K  |  ■  0  )  is  determined  by  plotting  |  K  |  vs.  X. 
For  the  test  problems  and  the  example  problems  discussed  herein  these  curves 
were  smooth  and  no  difficulties  were  encountered  in  implementing  this  scheme. 


I 


I 
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3.  TEST  PROBLEMS 


The  performance  of  the  8-node  element  was  tested  by  using  it  to  solve  two 
different  problems  with  known  solutions.  These  tests  and  their  results  are  de¬ 
scribed  below. 

3.1  3-D  CRACK  IN  ISOTROPIC  MATERIALS 

The  problem  of  stress  singularities  at  the  tip  of  a  quarter  infinite  crack  in 
an  isotropic  half  space  with  the  crack  front  normal  to  the  free  surface.  Fig. 
2a,  was  first  studied  by  Benthem  [10],  Using  a  semi  ana  1 yt i ca 1 -numer i ca 1  tech¬ 
nique  he  evaluated  the  eigenvalue  X  corresponding  to  a  state  of  symmetric  defor¬ 
mation  for  different  values  of  the  Poisson's  ratio  *.  These  results  were  subse¬ 
quently  verified  by  Bazant  and  Estenssoro  [9]  who  proposed  the  variational 
principle  used  in  the  present  study  and  developed  a  4-node  quadrilateral  element 
for  isotropic  materials.  In  addition  they  also  reported  the  eigenvalues  for  an¬ 
tisymmetric  deformations.  Later  Benthem  [11]  verified  these  results  using  a  fi¬ 
nite  difference  scheme.  Another  study  of  the  same  problem  has  been  reported  by 
Kawai  et  al.  [13].  While  confirming  the  results  of  [10]  they  also  reported  the 
detection  of  additional  stronger  singularities.  But  that  has  not  been  verified 
in  subsequent  work  [9]  and  now  it  is  generally  believed  that  the  results  of  Ben¬ 
them  [11]  and  Bazant  and  Estenssoro  [9]  are  correct  [3].  In  our  first  test  the 
8-node  element  was  used  to  solve  the  same  problem  and  the  results  were  compared 
against  those  from  [9,11]. 

The  domain  of  interest  in  the  [8 ,4>)  plane  is  a  rectangle  bounded  by  6m 0, 
0*ir/2,  <t>~0  and  Fig.  2b.  Only  one  half  of  the  problem  was  analysed  using 

the  symmetry  of  the  geometry.  The  domain  was  sub-divided  into  (n$  x  n ^  el¬ 
ements;  n^  in  the  8  direction  and  n^  in  the  <t>  direction.  At  nodes  on  the 
boundary  AB  corresponding  to  the  plane  of  symmetry  <£*ir,  the  symmetric  or  anti- 
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symmetric  boundary  conditions  were  imposed  depending  on  the  case  to  be  analysed. 
The  symmetric  boundary  conditions  are:  w=0,  crg^O,  or<£«0  and  the  anti -symme¬ 
tric  boundary  conditions  are:  u*0,  v=0,  <7^* 0.  Of  these  only  the  displace¬ 

ment  boundary  conditions  are  imposed  explicitly.  The  stress  free  conditions  are 
the  natural  boundary  conditions  of  the  variational  principle  and  will  be  auto¬ 
matically  satisfied  if  the  corresponding  displacement  component  is  not  pre¬ 
scribed.  No  displacement  boundary  conditions  are  prescribed  on  the  other  bound¬ 
aries  (i .e.  0A,  BC  and  CO  in  Fig.  2b)  so  that  the  corresponding  surfaces  remain 

traction  free. 

The  work  of  Bazant  and  Estenssoro  [9]  shows  that  the  4-node  quadrilateral  el¬ 
ement  based  on  the  present  variational  principle  would  produce  the  same  results 
as  those  reported  in  [10,11]  for  this  problem.  As  the  approaches  used  in 
[10,11]  are  quite  different  from  that  in  [9]  it  strongly  suggests  that  those  re¬ 
sults  are  accurate.  Therefore  it  is  reasonable  to  expect  that  the  8-node  ele¬ 
ment  too  would  yield  the  same  results.  The  more  important  aspect  of  this  test 
was  to  compare  the  errors  and  the  convergence  of  the  8-node  element  against 
those  for  the  4-node  element.  However,  the  first  task  was  to  decide  upon  the 
number  of  integration  stations  ng  to  be  used  for  this  problem  and  for  that 
purpose  one  particular  case  (corresponding  to  #'■■0.3  and  ant i -symmetr  i c  deforma¬ 
tion)  was  run  with  varying  values  for  ng..  The  results  are  presented  in  Table 
1.  On  the  basis  of  the  rapidity  of  convergence  to  a  steady  value  nQ*4  was 
chosen  as  the  optimal  order  of  integration  for  this  problem  and  was  used  in  the 
subsequent  runs. 

The  results  presented  in  Table  1  suggest  that  a  mesh  of  (3x6)  should  produce 
values  for  X  within  1%  of  the  'correct1  value.  Therefore  in  the  next  step  of 
computing  values  of  X  corresponding  to  different  values  of  v  the  mesh  of  (3x6) 
was  used.  As  an  additional  check  though,  a  mesh  of  (4x8)  also  was  used.  Subdi- 

-  19  - 


visions  were  made  uniform  since  it  is  reported  in  [9]  that  non  uniform  grids  are 
not  more  effective  than  uniform  grids.  The  number  of  elements  in  the  4>  direc¬ 
tion  (n,p  was  always  kept  equal  to  twice  that  in  the  8  direction  (n$)  so 
that  the  finite  elements  would  have  equal  side  lengths  on  the  (9,4>)  plane. 

For  the  above  two  meshes  and  different  values  of  v,  the  computed  values  of  X 
for  both  symmetric  and  ant i -symmetr i c  deformations  are  presented  in  Table  2. 
The  corresponding  'extrapolated  values'  from  [9]  and  the  values  from  [11]  are 
also  presented  in  Table  2  for  comparision  and  excellent  agreement  is  evident. 
What  is  more  important  though  is  the  remarkably  rapid  rate  of  convergence  exhib¬ 
ited  by  the  8-node  element  when  compared  with  that  for  the  L-node  element.  To 
see  this  the  results  in  Table  2  have  to  be  compared  with  those  in  Fig.  8  of  [9]. 
Comparision  should  be  based  on  the  total  number  of  degrees  of  freedom  (which  is 
the  same  as  the  total  number  of  equations)  because  that  is  the  primary  factor 
that  determines  the  computational  effort  involved.  As  a  result  of  this  rapid 
convergence  the  need  to  extrapolate  the  results  using  convergence  characteris¬ 
tics  reported  in  [9]  can  be  eliminated  and  reasonably  accurate  results  can  be 
directly  obtained  from  a  relatively  coarse  finite  element  mesh. 

3.2  AX  1 SYMMETR 1C  NOTCH  J_N  TRANSVERSELY  ISOTROPIC  MATERIAL 

The  second  test  problem  was  chosen  with  the  intention  of  studying  the  per¬ 
formance  of  the  element  (and  in  fact  of  the  whole  finite  element  scheme)  when 
the  material  is  anisotropic.  The  geometry  is  assumed  to  be  axi symmetric  and  the 
conical  notch  is  defined  by  6  ■  90,  Fig.  3*  The  region  0^9<>9  0  is  occupied 
by  the  material  which  is  transversely  isotropic  with  respect  to  the  z  axis.  In 
[6]  analytical  solutions  are  obtained  for  this  problem  for  both  traction  free 
and  rigidly  clamped  boundary  conditions  on  the  surface  9  -  9Q.  When  the  ma¬ 
terial  is  transversely  isotropic,  there  are  five  independent  elastic  constants 


for  Cjj.  For  a  test  problem  we  have  selected  the  following  special  class  of 
transversely  isotropic  materials  in  which  Cjj  depends  on  four  parameters  m.  >'» 

3  and  7: 

C„  -  (6  +  2m)  31 

C33  -  (5  +  2m)/  31 

C44  “ 

C,a  “  5 

C11  “  C12  " 

where 

5  -  2m*'/  (1  ■2«')  (i»3) 

This  class  of  materials  has  the  property  that  the  'eigenvalues  of  the  elasticity 
constants'  are  repeated  [15] •  Moreover,  when  y  *  3  “  1  we  have  isotropic  ma¬ 
terials.  We  chose  0Q  *  $ir/h  and  considered  the  traction-free  boundary  condi¬ 
tion. 

In  the  ( 6,<t> )  plane  the  materia)  occupies  a  rectangular  domain  given  by 

(0 0£<££2?r)  with  the  two  boundaries  <t>*0  and  <p*2tr  having  the  same  dis¬ 
placements.  In  general  the  finite  element  mesh  would  have  to  cover  this  area. 
However,  the  axi symmetry  of  the  present  problem  can  be  used  to  significantly  re¬ 
duce  the  magnitude  of  the  computational  effort.  On  any  plane  through  the  z  axis 
symmetric  boundary  conditions  (viz.  w-0,  og^o^-O)  should  be  satisfied. 
Further,  the  non-zero  displacements  u  and  v  should  be  independent  of  4>.  This 
enables  us  to  confine  the  analysis  to  just  one  strip  of  elements  in  the  8  direc¬ 
tion.  The  number  of  elements  in  that  strip  is  n$.  Through  appropriate  assem¬ 
bly  procedures  the  same  nodal  degree  of  freedom  is  assigned  to  the  value  of 
f{8,<t>)  or  g  {8  ,<t>)  at  all  nodes  at  the  same  4>  level.  Symmetric  boundary  condi¬ 
tions  are  prescribed  at  the  two  sides;  and  at  the  top  (i.e.  at  8*0)  the  axisym- 


metric  conditions  (viz.  u  “0,  w*0)  are  imposed.  The  bottom  boundary  0*9o  is 
left  traction  free.  The  width  of  the  strip  of  elements  is  adjusted  according  to 
n$  so  as  to  result  in  square  elements. 

The  material  constants  C j j  can  be  normalized  by  dividing  by  m.  Therefore 
the  results  would  depend  only  on  v and  y.  In  the  first  case,  for  a  particular 
set  of  their  values  (f“0.3,  3*1 -5,  V*0.5)  X  was  computed  using  different  numbers 
of  integration  stations  and  different  meshes.  The  results  are  presented  in  Ta¬ 
ble  3>  As  described  in  the  case  of  test  problem  1,  in  this  case  too  an  optimal 
value  for  nG  was  selected  as  4.  A  mesh  with  4  elements  gives  results  accurate 
within  0.1%.  This  combination  was  then  used  to  compute  values  of  X  correspond¬ 
ing  to  other  values  of  v,  3  and  y  and  the  results  are  given  in  Table  4  along 
with  the  analytical  solutions  from  [6].  The  special  case  of  0  ■  7  ■  1  corre¬ 
sponds  to  isotropic  materials  and  the  value  of  X  for  that  case  has  been  reported 
in  [16].  The  agreement  between  the  finite  element  results  and  the  analytical 
results  is  quite  encou-aging  and  reinforces  the  belief  that  the  finite  element 
scheme  and  the  8-node  element  would  perform  satisfactorily  even  in  the  case  of 
anisotropic  materials.  Once  again,  the  remarkable  rapidity  with  which  the  re¬ 
sults  converge  as  evident  from  Table  3  should  be  noted. 


4.  AN  APPLICATION  TO  COMPOSITE  MATERIALS 
The  use  of  light  .weight  composite  materials  in  industrial  applications  is 
steadily  increasing.  In  areas  where  weight  is  a  crucial  factor  such  as  space¬ 
craft  designs,  one  would  like  to  be  able  to  exploit  the  strength  of  the  materi¬ 
als  to  the  fullest  possible  extent  and  the  availablity  of  relevant  design  infor¬ 
mation  assumes  great  importance.  Failure  criteria  are  among  those  highly 
desired  material  parameters.  In  the  case  of  laminated  composites  experimental 
observations  indicate  that  failure  may  occur  along  the  interface  between  the 
layers  or  transverse  to  the  layers.  It  is  well  known  that  stress  singularities 
occur  at  the  free-edge  of  the  interface  between  two  dissimilar  materials  and 
they,  no  doubt,  play  an  active  role  in  initiating  these  failures.  As  such,  de¬ 
velopment  of  failure  criteria  requires  a  sound  understanding  of  the  nature  of 
stress  singularities  that  would  be  present  at  the  interface  and  it  has  been  the 
subject  of  many  recent  investigations,  (e.g.  see  [  1 7“ 1 9] ) • 

Consider  the  case  where  a  transverse  crack  is  already  present  at  the  free 
surface  of  a  laminated  composite.  It  may  have  resulted  from  a  material  imper¬ 
fection  or  from  a  fabrication  error.  Or  else  it  may  have  been  initiated  by  the 
above  mentioned  free-edge  stress  singularities  arising  as  a  result  of  external 
loading.  In  Fig.  4  such  a  crack  which  is  normal  to  and  ends  at  the  interfaces 
is  shown  by  the  shaded  area.  The  presence  of  this  crack  would  affect  the  local 
stress  field  and  the  nature  of  stress  singularities.  In  fact  it  would  give  rise 
to  a  new  stress  singularity  (i.e.  in  addition  to  those  that  would  occur  along 
free-edges  such  as  MN)  along  the  transverse  crack  edge  MQ.  A  similar  situation 
would  occur  with  any  other  orientation  of  the  crack,  and  in  particular  with  a 
crack  along  the  interface  which  could  result  from  delamination.  But  for  the 
purpose  of  the  present  discussion  we  shall  confine  our  attention  to  the  crack 
shown  in  Fig.  4.  The  singularity  along  MQ  would  now  control,  at  least  locally, 
the  process  of  material  failure  and  as  such  merits  investigation. 
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Even  though  in  principle  the  singularity  at  a  point  on  MQ  is  3"dimensiona1 , 
as  a  first  approximation  its  analysis  may  be  reduced  to  a  2-  dimensional  problem 
provided  that  the  point  is  sufficiently  away  from  the  free  surface.  This  has 
been  done  in  [19]  and  also  given  therein  is  an  account  of  relevant  previous  work 
on  2-dimensional  problems.  In  the  present  study  we  direct  our  attention  to  the 
nature  of  the  3"din>ensional  stress  singularity  at  the  point  M. 

In  the  immediate  neighbourhood  of  the  point  M  the  problem  of  stress  singular¬ 
ities  is  very  much  similar  to  the  problem  of  quarter  infinite  crack  in  a  half 
space  discussed  in  Section  3  except  for  the  fact  that  the  material  is  neither 
homogeneous  nor  isotropic.  The  region  (x<0,  z£0)  is  occupied  by  the  material  of 
layer  2  while  the  regions  (x>0,y>0,z&0)  and  (x>0,y<0,z£0)  are  occupied  by  the 
material  of  layer  1.  It  should  also  be  noted  that  even  though  in  Section  3  it 
was  possible  to  confine  the  analysis  to  only  one  half  of  the  problem  by  exploit¬ 
ing  symmetry,  in  the  present  case  the  material  will  not  exhibit  such  symmetry  in 
general  and  therefore  the  full  problem  has  to  be  analysed.  '  Once  these  points 
are  identified  the  analysis  can  be  performed  in  a  mannar  quite  parallel  to  that 
employed  previously. 

The  domain  occupied  by  the  material  in  (9  A)  plane  is  the  rectangle  (0£9£ir/2, 
0«t><2n)  .  In  this,  the  two  regions  (0£8$ir/2,  0«t><ir/2)  and  (O£0£?r/2,  }rr/2«i><2ir) 

correspond  to  material  1  and  the  region  (OS0S7r/2,  n/2<4><}ir/2)  corresponds  to  ma¬ 
terial  2.  The  entire  region  is  subdivided  into  (ng  x  n,p  finite  elements. 

For  this  problem  n^  was  made  equal  to  4ng  so  that  square  elements  would  re¬ 
sult.  It  also  ensures  that  each  element  would  lie  entirely  within  one  material. 
During  the  computation  care  has  been  exercised  in  assigning  the  correct  material 
properties  to  each  finite  element. 

The  remaining  steps  are  routine.  The  element  stiffnesses  are  computed  for  a 
given  X  and  assembled  into  the  global  stiffness  matrix.  The  boundary  conditions 
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are  all  traction  free  conditions  and  none  needs  to  be  imposed  explicitly.  The 
eigenvalue  X  is  obtained  by  plotting  the  determinant  of  the  global  stiffness  ma¬ 
trix  against  corresponding  value  of  X. 

As  explained  in  Section  3  the  first  task  is  to  select  the  optimal  order  of 
Gauss  integration  ng.  Like  in  the  test  problems,  this  is  done  by  carrying  out 
numerical  tests  to  study  the  convergence  of  the  results  as  the  mesh  is  refined 
and  ng  is  varied.  Once  optimal  ng  is  selected  a  mesh  that  could  be  expected 
to  yield  'sufficiently  accurate1  results  can  be  selected  and  used  for  analysing 
the  remaining  cases  of  the  same  type.  What  is  'sufficiently  accurate'  is,  of 
course,  subjective  and  would  depend  to  a  great  extent  on  the  intended  applica¬ 
tion  of  the  result.  It  would  also  be  controlled  by  physical  limitations  of 
available  computing  capacity. 

For  the  purpose  of  a  numerical  example,  we  assume  that  each  layer  in  the  com¬ 
posite  is  made  of  a  fiber  reinforced  material  which  can  be  regarded  as  an  ortho- 
tropic  material  whose  material  symmetry  is  with  respect  to  the  x,,  x2,  x3 
axes  where  x^x  and  x3  axis  is  the  fiber  direction  which  makes  an  angle  a 
with  the  negative  y  axis  as  shown  in  Fig.  4.  We  further  assume  that  each  layer 
is  made  of  the  same  material  although  the  ply  angle  a  may  vary  from  layer  to 
layer.  The  value  of  a  corresponding  to  layer  1  is  a(  and  that  corresponding 
to  layer  2  is  a2.  Referring  to  the  (xlfx2,x3)  axes  let  the  orthotropic 
material  have  the  following  material  properties  correspond i ng  to  T300/5208  gra¬ 
phite/epoxy  given  in  [20].  This  is  the  same  as  the  material  referred  to  as  com¬ 
pos  i te  T  in  [18,19]. 
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are  the  shear  moduli 


In  the  above  equation  E;  are  the  Young's  moduli,  G;j 
and  Vjj  are  the  Poisson's  ratios  [21].  They  can  be  used  to  compute  the  ma¬ 
terial  compliance  matrix  [21,22]  which  when  inverted  using  the  relations  given 
in  [21]  yields  the  material  stiffness  matrix  referred  to  the  (x,,x2,x  ) 
axes.  Then  the  material  stiffness  matrix  C j j  referred  to  the  (x,y,z)  axes  can 
be  obtained  by  an  appropriate  co-ordinate  transformation.  The  use  of  ply  angle 

ax  will  result  in  Cjj  corresponding  to  material  1  while  a}  will  result  in 
Cjj  for  material  2. 

For  the  purpose  of  selecting  an  optimal  value  for  nQ  we  used  the  particular 
case  with  c^-ir/2,  a2*0.  This  choice  restored  material  symmetry  with  respect 
to  the  y*0  plane  and  thus  allowed  the  analysis  of  only  one  half  of  the  problem 
as  in  Section  3-  That  way  the  refinement  of  the  mesh  could  be  carried  further 
than  it  would  have  been  possible  if  the  full  problem  had  to  be  analysed.  Atten¬ 
tion  was  focused  only  on  the  smallest  admissible  value  of  X.  The  results  are 
presented  in  Table  5.  On  the  basis  of  rapid  convergence  to  a  steady  value 

■  6  was  chosen  as  being  optimal  for  this  type  of  problem.  The  need  for 
more  integration  stations  than  that  required  in  the  isotropic  case  discussed  in 
Section  3  is  perhaps  a  reflection  of  the  higher  degree  of  anisotropy  involved. 
However,  examination  of  the  results  in  Table  5  indicate  that  even  in  this  highly 
anisotropic  problem  the  numerical  scheme  performs  satisfactorily. 

A  mesh  of  (3x12)  was  selected  for  use  in  ar^lysing  other  (ay/a2)  combina¬ 
tions.  They  would  lack  the  material  symmetry  about  the  plane  y»0  required  to 
confine  the  analysis  to  one  half  of  the  problem  and  would  necessitate  the  analy¬ 
sis  of  the  full  problem.  The  available  computing  capacity  did  not  permit  the 
use  of  a  more  refined  mesh.  However,  based  on  the  results  of  the  above  test 
problem  and  also  on  the  convergence  pattern  displayed  by  the  results  correspond¬ 
ing  to  other  ply  angle  combinations  we  feel  that  the  results  from  the  (3x12) 


mesh  would  be  accurate  to  within  at  least  2%.  In  fact  in  most  of  the  cases  the 
error  would  be  less  than  1%.  Proper  verification  of  that  estimate  requires  a 
reliable  analytical  result  or  at  least  a  result  from  a  different  numerical  pro¬ 
cedure  either  of  which  to  the  authors'  knowledge  is  not  available  as  yet.  The 
results  from  a  mesh  of  (3x12)  for  different  ply  angle  combinations  are  presented 
in  Table  6.  For  comparision,  in  the  same  table  we  have  provided  within  paren¬ 
thesis  the  results  from  a  mesh  of  (2x8).  The  small  differences  between  the  two 
sets  of  results  lend  support  to  the  belief  that  the  results  are  reasonably  accu¬ 
rate.  In  view  of  the  symmetry  property  that  the  X  corresponding  to  (a)/a2) 
is  the  same  as  that  for  (-a^/-a2)  only  the  values  of  X  for  OSaiSjr/2  are  pre¬ 
sented.  They  vary  from  0.3016  (for  0/90  combination)  to  0.4572  (for  90/60  com¬ 
bination).  The  strongest  singularity  occurs  for  the  combination  (0/90)  and  the 
corresponding  order  of  stress  singularity,  1-X,  is  -0.6984.  In  [19]  the  two-di¬ 
mensional  singularity  at  a  point  sufficiently  away  from  the  free  surface  was 
studied  and  it  is  shown  that  the  strongest  singularity  in  that  case  also  occurs 
for  the  (0/90)  combination. 
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.  STRESS  SINGULARITIES  AT  OTHER  THREE-DIMENSIONAL  CORNERS 


As  a  further  example,  in  this  section  the  present  numerical  scheme  is  used  to 
investigate  the  possible  occurence  of  stress  singularities  in  other  3“D  corners. 
The  problem  of  a  corner  formed  by  three  mutually  perpendicular  planes  is  consid¬ 
ered  in  the  computations.  The  three  planes  are  defined  by  <t>  ■  0;  <t>  ■  7r/2  and 

6  *  jt/2.  (See  Figure  5).  Two  different  situations  arise  depending  on  which 
part  is  occupied  by  material  and  which  part  is  void.  Let  the  region 
(  0  £  <t>  £  ir/2,  0  <  9  £  jt/2  )  be  denoted  by  R.  One  problem  is  when  R  is  void  and 
the  rest  of  the  space  is  occupied  by  the  material.  Hereafter  it  will  be  refer¬ 
red  to  as  the  'exterior  problem'.  The  other  problem  arises  when  the  situation 
is  reversed,  i.e.,  R  is  occupied  by  the  material  while  the  rest  is  void.  It 
will  be  referred  to  as  the  'interior  problem'. 

Two  materials  are  used  in  the  sample  computations:  isotropic  material  (with 
Poisson's  ratios  of  0  and  0.3)  and  the  material  referred  to  as  composite-T  in 
Section  4  with  properties  given  by  equation  (44).  In  the  latter  case  the  x3 
axis  is  taken  to  be  in  the  z  direction  so  that  material  properties  are  rotation- 
ally  symmetric  about  the  z  axis.  This  enables  the  analysis  to  be  confined  to 
one  half  of  the  problem  and  be  performed  separately  for  symmetric  and  antisymme¬ 
tric  modes  of  deformation. 

5.1  EXTERIOR  PROBLEM 

One  half  of  the  domain  of  interest  is  shown  in  Figure  6(a).  The  number  of 
elements  in  the  9  direction  between  0  and  w/2  is  the  same  as  that  between  n/2 
and  ir  and  is  denoted  by  ng.  The  number  of  elements  in  the  direction  between 
0  and  ir/4  is  n,^  while  that  between  7r/4  and  ir  is  n^.  In  order  to  have 

square  e  «  ients  on  the  9-<t>  plane  the  ratio  ngtn^sn^  was  maintained  at 
2:1:3- 


-•  -• j  *  ,  t  *-  ^  -  *- ■ '  -  -  *  -  -  -  - 


The  coarse  mesh  had  (ng,  n^ ,  n^)  ■  (2,1,3)  and  the  refined  mesh 

(14,2,6).  Based  on  the  experience  gained  from  the  computations  of  the  previous 
problems  it  is  believed  that  the  latter  mesh  would  produce  results  within  2%  of 
the  exact  value.  However  the  large  size  of  the  domain  prevented  any  further  re¬ 
finements.  The  coarse  mesh  was  first  used  to  scan  the  real  axis  from  0  to  1.0 
in  search  of  a  root  X.  When  one  was  located  the  refined  mesh  was  used  to  com¬ 
pute  its  value  more  accurately.  Previous  experience  suggested  that  ng  -  5 
should  be  adequate  (where  ng  is  the  number  of  Gauss  integration  stations). 
But  for  purposes  of  comparision  calculations  were  performed  using  other  values 
of  ng  also.  The  results  obtained  for  both  isotropic  and  composite-T  materials 
are  presented  in  Table  7«  In  each  case  two  symmetric  and  one  ant i -symmetr i c 
roots  were  observed  on  the  real  axis  between  0  and  1.0.  They  all  would  give 
rise  to  stress  singularities. 

5.2  INTERIOR  PROBLEM 

Shown  in  Figure  6(b)  is  one  half  of  the  domain  of  interest  for  this  problem 
on  the  6-<t>  plane.  The  number  of  elements  in  the  d  direction  is  ng  and  that  in 
the  4>  direction  is  n^,.  Ratio  ng:n^,  was  maintained  at  2:1  so  that  square 
elements  would  result.  In  order  to  obtain  the  same  size  of  elements  as  in  the 

exterior  problem  (n^,n,p  *»  (2,1)  was  used  for  the  coarse  mesh  while  (L,2) 
was  planned  for  the  refined  mesh.  However  in  the  cases  of  both  isotropic  and 
composite-T  materials  a  search  using  the  coarse  mesh  revealed  that  there  are  no 
real  roots  of  X  capable  of  causing  stress  singularities  (i.e.  between  0  and 
1  .0)  . 

Attention  was  then  focused  on  complex  roots  of  X  to  see  whether  there  exists 
a  X  whose  real  part  is  between  0  and  1.0.  Computations  were  carried  out  using 
complex  arithmatic  and  appropriate  changes  of  the  relevant  variables  from  real 
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to  complex.  A  rectangular  region  on  the  complex  X  plane  was  scanned  for  roots. 
The  search  was  performed  by  subdividing  the  region  into  a  grid  of  small  rectan¬ 
gles  and  then  checking  each  of  them  to  see  whether  a  line  corresponding  to 
Re(|K|)  ■  0  and/or  lm(|K|)  *  0  would  be  present  inside  it.  (Re  and  Im  stand 
for  the  Real  part  and  the  Imaginary  part  respectively).  This  was  determined  by 
comparing  the  signs  of  the  relevant  quantity  (i  .e.  either  Re  ( | K I )  or  lm(|Kl)  as 
the  case  may  be)  at  the  four  corners.  If  the  sign  changes  among  the  four  cor¬ 
ners  then  it  would  be  concluded  that  the  corresponding  line  passes  through  the 
rectangle.  By  searching  each  small  rectangle  this  way,  it  was  possible  to  trace 
the  locus  of  the  lines  Re(|K|)"0  and  Im(|K|)*0.  Where  they  intersect  is  a  root 
of  X.  To  illustrate  the  type  of  plot  one  would  obtain,  the  results  correspond¬ 
ing  to  composite-T  material  are  given  in  Figures  7  and  8.  Within  the  region 
scanned  no  roots  capable  of  causing  stress  singularities  were  found  even  though 
some  roots,  both  real  and  complex,  were  present  in  the  region  Re  (X)  >  1.0. 
(Note  that  X  «  0  corresponds  to  a  rigid  body  translation  while  X  »  1  produces  a 
rigid  body  rotation  and  both  of  them  produce  no  stresses).  In  the  case  of  iso¬ 
tropic  materials  with  v»0  and  0.3  a  similar  situation  prevailed.  Correspond¬ 
ing  plots  of  Re(|K|)«0  and  lm(|K|)«0  were  qualitatively  similar  to  Figs.  7  and  8 
and  hence  are  not  presented  here.  They  can  be  found  in  [23]. 


6.  CONCLUDING  REMARKS 


Assuming  that  the  displacements  in  an  anisotropic  material  is  propor¬ 
tional  to  r^  a  finite  element  scheme  is  proposed  to  determine  the  eigenvalue 
X.  Since  the  stress  is  proportional  to  kr^”'  where  k  is  the  proportionality 
factor,  the  stress  is  singular  when  Re  (X)  <1.  X  is  a  root  of  the  determinant 

of  the  matrix  K  in  Eq.  (40).  The  proportionality  factor  k  has  to  be  deter¬ 
mined  by  considering  the  complete  boundary  conditions  which  include  the  boundary 
S,  in  Fig.  1.  If  k  happens  to  be  zero,  there  is  no  stress  singularity  at 

r  *  0  even  if  Re  (X)  <  1.  Thus  the  existence  of  a  singularity  is  not  certain  un¬ 
til  the  complete  boundary  conditions  are  considered.  If  a  singularity  exists, 
the  solution  obtained  here  provides  the  order  of  stress  singularity  and,  if  one 
desires,  the  stress  distribution  near  r  =0  except  the  proportionality  factor 
k . 

Complications  arise  when  X  is  a  multiple  root  of  the  determinant  of  K. 

*v 

The  related  problem  for  two  dimensional  cases  has  been  investigated  in  [24,25] 
for  isotropic  materials  and  in  [26]  for  anisotropic  materials.  It  was  shown 
that  in  addition  to  the  r^“^  singularity  the  stress  may  have  the 
r^'1  (In  r)  singularity.  The  conditions  for  which  the  stress  has  the 
r^-1  (In  r)  singularity  were  given  in  [24].  The  same  conditions  apply  to  the 
matrix  K  even  though  our  problem  is  a  three-dimensional  one.  No  multiple 

roots  are  found  in  the  examples  presented  here  nor  in  other  three-dimensional 
problems  reported  in  the  literature  except  in  [6],  Even  in  the  case  of  multiple 
roots  of  | K | ,  if  the  curves  Re(|K|)«0  and  lm(|K|)=0  are  plotted  on  the  complex 
plane  as  in  Figs.  7  and  8,  a  root  of  multiplicity  m  would  present  itself  as  a 
common  point  of  intersection  of  m  curves  of  Re(|K|)“0  and  m  curves  of 
lm(|K|)-0. 


Another  complication  arises  when  the  conical  wedge  surface  S,  in  Fig.  1, 
is  not  traction  free  and  the  matrix  K  is  singular.  In  this  case  the  right 

hand  side  of  Eq.  (hO)  is  non-zero  and  a  solution  may  not  exist.  A  modified  so¬ 

lution  would  give  a  singularity  of  the  form  k *(  1  n  r)  in  which  k  is  a  con¬ 
stant.  The  related  problem  for  two-dimensional  cases  has  been  studied  in 
[27-29]  for  isotropic  materials  and  in  [30]  for  anisotropic  materials.  It 
should  be  pointed  out  that,  unlike  the  constant  k  mentioned  earlier  which  is 
indeterminate  until  the  complete  boundary  conditions  are  considered,  k*  here 
can  be  determined  without  solving  the  complete  boundary  value  problem  [30] . 
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Appendix  A 


LIST  OF  COMPUTER  PROGRAMS  DEVELOPED  FOR  THIS  STUDY 

(1)  Program:  CON. NOTCH 

This  program  is  for  analysing  the  possible  stress  singularities  at  the 
apex  of  conical  regions  in  transversely  isotropic  materials  using  8-node 
quadrilateral  finite  elements.  The  conical  region  can  be  defined  either  as 
having  one  conical  boundary  or  as  having  two  conical  boundaries  with  ma¬ 
terial  in  between  them.  Either  traction  free  or  rigidly  clamped  conditions 
can  be  imposed  on  each  boundary.  The  user  has  control  over  the  number  of 
elements  and  the  order  of  integration.  The  program  can  be  used  either  to 
locate  a  real  root  in  a  specified  interval  on  the  real  axis  or  to  compute 
the  determinants  corresponding  to  given  real  values  of  X. 

(2)  Program; _ D3. CRACK. HALF 

This  is  for  the  analysis  of  possible  stress  singularity  at  the  tip  of 
the  transverse  crack  in  a  composite  laminate  shown  in  Fig.  L.  Material 
properties  are  assumed  to  be  such  that  there  is  symmetry  and  only  one  half 
of  the  problem  is  analysed.  8-node  quadrilateral  finite  elements  are  used. 
The  user  can  vary  the  mesh  details  and  the  order  of  Gauss  integration.  The 
program  is  capable  of  either  searching  for  a  real  root  within  a  given  in¬ 
terval  of  the  real  axis  to  a  specified  accuracy  or  of  just  computing  the 
determinants  for  specified  real  values  of  X. 


The  same  program  can  be  used  to  analyse  the  3-D  crack  in  an  isotropic 
half  space  also  by  assigning  isotropic  material  properties  to  both  layers 


of  the  composite. 

(3)  Program: _ D3. CRACK. FULL 

This  is  for  the  analysis  of  possible  3"D  stress  singularities  at  the  tip 
of  a  crack  in  a  laminated  composite.  No  material  property  symmetries  are 
assumed  and  the  full  problem  is  analysed.  The  crack  can  be  either  normal 
to  the  interface  or  along  the  interface.  Mesh  is  generated  using  8-node 
elements.  Mesh  details  and  the  order  of  Gauss  integration  are  user  defina¬ 
ble.  The  program  can  either  search  for  a  real  root  in  a  given  interval  of 
the  real  axis  to  a  specified  accuracy  or  just  compute  the  determinants  cor¬ 
responding  to  given  real  values  of  X. 

00  Program:  D3. CORNER. EXT 

The  exterior  problem  for  a  3“D  corner  formed  by  3  mutually  perpendicular 
planes  is  analysed  for  the  possible  occurence  of  stress  singularities  by 
this  program.  Mesh  is  generated  for  8-node  quadrilateral  elements.  Mesh 
parameters  and  the  order  of  Gauss  integration  are  user  definable.  Material 
is  required  to  have  properties  symmetric  with  respect  to  the  plane  of  geo¬ 
metric  symmetry  for  the  problem.  Analysis  is  performed  for  either  the  sym¬ 
metric  or  ant i -symmetr i c  mode  and  the  program  can  either  search  for  a  real 
root  within  a  specified  region  on  the  real  axis  to  a  specified  accuracy  or 
compute  the  determinants  corresponding  to  given  real  values  of  X. 

(5)  Program: D3 . CORNER . I  NT 

Similar  to  D3-C0RNER.EXT  except  that  this  analyses  the  interior  problem. 
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(6)  Program: _ D3. CORNER. INT.CMPLX 

Similar  to  D3-C0RNER. I  NT  except  that  the  search  for  a  root  is  carried 
out  on  the  complex  plane.  Program  produces  the  plots  of  the  locii  of 
Re(|K|)  ■  0  and  of  lm(|K|)  -  0  in  a  specified  region  on  the  complex 
plane.  Points  of  intersection  of  the  two  types  of  lines  indicate  locations 
of  roots.  An  option  to  refine  the  root  upto  a  specified  accuracy  also  is 
avai table. 


Number  of 
Gauss 

Integration 


Mesh  ng  x  n^ 
(Total  number  of  DOF) 


Stations 

2  x  4 

3  x  6 

4x8 

nG 

(111) 

(219) 

(363) 

3 

0.3890 

O.3883 

0.3879 

4 

O.39I.6 

0.3912 

0.3911 

5 

0.3970 

0.3929 

O.3920 

6 

0.3991 

0.3941 

0-3932 

8 

0.4014 

0.3955 

0.3943 

10 

0.4029 

0.3964 

0.3950 

*Mesh  is  foi 

one  half  of  tl 

ie  problem. 

[Mllll 


value 

of 


Symmetric  Deformation 


Anti -symmetric  Deformation 


8-node  < 
nc  -  A 

si ement 

from 
Bazant  & 
Estenssoro 

from 

Benthem 

cm 

8 -node 
n6  ■  4 

element 

mesh 

(3x6) 

mesh 

(4x8) 

[9] 

mesh 

(3x6) 

mesh 

(4x8) 

0.50J0 

0.5003 

0.5 

0.5 

O.4985 

0.4983 

0.5107 

0.5099 

- 

- 

0.4492 

0.4492 

0.5173 

0.5166 

0.5164 

0.5164 

0.4316 

0.4315 

0.5487 

0.5478 

0.5477 

0.5477 

0.3912 

0.39H 

0.5883 

O.5869 

0.5868 

0.5868 

■ 

0.3703 

0.3701 

®  |  from 

t  &  I  Ben them 

cm 


Number  of 
Gauss 

Integration 

Stations  n_ 

G 

Number  of 

Elements  ng 

2 

3 

4 

5 

6 

3 

0.5999 

0.5309 

0.5292 

0.5299 

0.5298 

u 

0.5881 

0.5306 

0.5293 

0.5299 

0.5298 

5 

0.5880 

0.5307 

0.5293 

0.5299 

0.5298 

Table  3  Values  of  X  for  a  conical  notch  (Fig.  3,  0o*3n/4) 

in  a  tranversely  isotropic  material  (v*0.3,  6*1.5, 
y*l,  y*0.5)  using  8-node  element  with  different 
numbers  of  Gauss  integration  stations  and  different 
meshes.  Exact  solution  is  X* 0.5296  [6]. 
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Material  Properties 

Using  8-node 
finite  elements 

nG-4 

Analytical 
solution 
from  [6] 

V 

3 

— 

y 

0.3 

1  .0 

1 .0 

0.8014 

0.8012* 

1  .0 

2.0 

0.9445 

0.9437 

1.5 

0.5 

0-5293 

0.5296 

1.5 

2.0 

0.8529 

0.8533 

0.5 

0.5 

0.9160 

0.9131 

*this  case  corresponds  to  isotropic  material  and 
agrees  with  the  analytical  solution  for  isotropic 
cones  reported  in  [ 1 6] 


Table  4  Comparison  of  X  for  the  conical  notch.  Fig.  3,  in 
a  transversely  isotropic  material  using  8-node  finit 


Humber  of 
Gauss 

Mesh*  ng  > 

integration 
Sea Cions 

2x8 

3x12 

4x16 

5x20 

3 

0.4220 

0.4353 

0.4393 

0.4412 

4 

0.4394 

0.4426 

0.4448 

0.4458 

5 

0.4448 

0.4469 

0.4482 

0.4484 

6 

0.4487 

0.4498 

0.4505 

0.4502 

8 

0.4520 

0.4533 

0.4533 

0.4524 

10 

0.4555 

0.4555 

0.4550 

0.4538 

16 

0.4597 

0.4589 

0.4578 

0.4560 

24 

0.4623 

0.4611 

0.4595 

0.4575 

0 . 3860 
(0.3909) 

0.35U 

(0.3557) 

0.3171 

(0.3181) 

0.3016 

(0.2998) 


o 

o 

rr\ 

O 

O 

sO 

90° 

0.4048 

0.4416 

(0.4089) 

(0.4494) 

0.3881 

0.4352 

(0.3904) 

(0.4377) 

0.3872 

0.4128 

0.41*97 

(0.3914) 

(0.4190) 

(0.4484) 

0.39H 

0.4020 

0.41*77 

(0.3974) 

(0.409D 

(0.41*63) 

0.4067 

0.4264 

0.4572 

(0.4128) 

(0.4356) 

(0.4602) 

0.4115 

0.4293 

0.4550 

(0.4206) 

(0.4362) 

(0.451*9) 

Table  6  Values  of  smallest  X  for  [ay/a2)  composite. 

Fig.  4,  obtained  from  a  (3x12)  mesh.  X  obtained 
from  a  (2x8)  mesh  are  shown  in  parentheses. 


Mode 


Mater ial 


Number  of 
Gauss 

Integration 

Stations 


1 sotropic 

\)  -  0 

1 sotropic 
v  *  0.3 

Compos i te 
-T 

0.7118 

(0.7254) 

0.7660 

(0.7803) 

0.7312 

(0.7457) 

0.7122 

(0.7269) 

0.7668 

(0.7824) 

0.7315 

(0.7440) 

0.8554 

(0.8643) 

0.7805 

(0.7928) 

0.8559 

(0.8654) 

0.7809 

(0.7934) 

0.8165 

(0.8356) 

0.7140 

(0.7294) 

O.7667 

(O.7807) 

0.7623 

(0.7680) 

0.7141 

(0.7301) 

0.7671 

(0.7823) 

0.7663 

(0.7744) 

Symmetr i c 


Symmetr i c 


Ant i -symm. 


10 


10 


10 


Table  7  Values  of  real  X  within  (0,  1.0)  for  3-D  corner 
exterior  problem.  Fig.  5a,  computed  on  a  (4,2,6) 
mesh.  Results  from  a  (2,1,3)  mesh  are  shown  in 
parentheses. 


3-D  body 


cross  section 


Fig.  3  Ax i symmetric  conical  notch 
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SYMMETRIC 


e(|K|)  *  0  and  lm(|K|)  ■  0  on  complex  X  plane  corresponding 
corner  interior  problem,  Fig.  5b,  for  composite  -  T  material 
i-symmetric  mode  of  deformation. 
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